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Abstract 



The relevance of various residue positions for the stability and the folding 
characteristics of the prion protein are investigated by using molecular dynam- 
ics simulations of models exploiting the topology of the native state. Highly 
significant correlations are found between the most relevant sites in our anal- 
ysis and the single point mutations known to be associated with the arousal 
of the genetic forms of prion disease (caused by the conformational change 
from the cellular to the scrapie isoform). Considerable insight into the con- 
formational change is provided by comparing the folding process of prion and 
doppel (a newly discovered protein) sharing very similar native state topol- 
ogy: the folding pathways of the former can be grouped in two main classes 
according to which tertiary structure contacts are formed first enroute to the 
native state. For the latter a single class of pathways leads to the native state. 
Our results are consistent and supportive of the recent experimental findings 
that doppel lacks the scrapie isoform and that such remarkably different be- 
havior results from differences in the region containing the two (5— strands and 
the intervening helix. 



I. INTRODUCTION 

Neurodegenerative diseases causing transmissible spongiform encephalopathies (TSE) 
are the subject of intense research. They affect humans and animals and include scrapie 
in sheep, mad-cow disease in bovines and Creutzfeldt- Jacob disease in humans. Most of 
them arise sporadically, the rest of the instances are inherited or transmitted by inoculation 
through the dietary or infected tissues. 

The belief is that the normal, and benign, form of the prion protein PrP c , which is found 
in various organs of vertebrates, including the brain, can undergo a post translational process 
leading to a conformational change of its native state |]|. This new form is designated PrP sc 
and contrarily to PrP c it is insoluble. The NMR investigation of the three-dimensional 
structure of recombinant PrP c of various species pHTOf has revealed that the N-terminus is 



unstructured whereas the C-terminus consists of three a-helices and two short /3-strands. 

On the other hand the structure of the malign form, PrP sc , is mostly unknown. However, 
spectroscopic investigations, mainly based on circular dichroism, showed that PrP c contains 
about 40% a-helices and little /3-sheet, whereasPrP sc is composed of about 30% a-helices 



and 45% /3-sheet [y,y 



Unfolding and re-folding experiments suggest that the conversion between the a-rich 
and /3-rich form occurs through a complete unfolding of the protein [ 13| . Experiments and 



clinical tests suggest two main processes for the onset and spreading of the prion infection. 
In the first case, the spreading of the scrapie form in healthy tissues contaminated with 
p r psc revea ] s that the presence of PrP sc helps further conversion from normal to scrapie 
form, i.e. PrP sc acts as a template for the restructuration of PrP c . In the second case, the 
genetic influence on the propensity of cellular PrP c to rearrange into PrP sc is revealed by 
the discovery of, at least, 20 single point mutations in the PrP gene in Humans which favor 
the spontaneous onset of the disease |]4],|5|JT4]-[T6l. 
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Within the picture presented in |13[] the role of the mutations in the arousal of the disease 
could be caused by several factors ranging from the lowering of a free energy barrier in the 
conversion process to the increase in the oligomerization rate. One of the goals of the present 
study is to understand which are the sites where mutations are expected to have the major 
impact in the stabilization/destabilization of the PrP c structure and/or its folding process. 

Recently, some experiments on transgenic mice lacking the PrP gene (PrnP Q /°) have 
revealed the onset of neurodegeneration (Purkinje-cell death) [T7J completely different from 
usual prion disease. This disease was ultimately traced back to a prion homologous gene, 
Prnd, (not expressed in normal subjects), which encoded a protein named doppel, Dpi 
which has almost the same topology of PrP c , despite the low sequence identity (25%). 

A striking difference between Dpi and PrP sc , is that the former affects the central 
nervous system while retaining its native structure, i.e. without a need to convert to a 



scrapie-like conformation. Various hypothesis have been formulated |L8[ to explain these 
different behaviours. 

The theoretical framework we elaborated and adopted in our study helps to shed light 
on the putative causes for such remarkably different behavior of doppel and prion. Indeed, 
the different folding behavior and stability of prion and doppel are related, at least in part, 
to the small, but important, differences in the topology of the native state which affect in an 
amplified fashion, the folding routes and ultimately facilitating structural rearrangements of 
prion. 

The model we adopt builds on the importance of the native state topology in steering 
the folding process, that is in bringing into contact pairs of amino acids that are found 
in interaction in the native state. In the past few years an increasing amount of experi- 
mental [|P9H 2"T| and theoretical |22|-30] evidence (such as all-atom molecular dynamics (MD) 
simulations in implicit solvent) have confirmed this view pi],P2 • 



A further confirmation of this view is provided by the remarkable accord of the key folding 
stages predicted by topology-based models and the available experimental results PB| , |3"1 



Indeed, the success of the topological picture itself helps to explain why the folding process 
is not too sensitive to the detailed chemical composition of most residues in a protein. The 
dependence of the folding process on the detailed chemistry is much more subtle than that 
given by the native state topology. 

Usually the folding mechanism is affected only when mutations occurs in a small set of 
key residues. Within our theoretical framework those key residues take part to contacts that 
are crucial for the folding process |I3 . 



The establishment of the key contacts leads to a rapid formation of further interactions. 
Interestingly, the folding bottlenecks can be identified by just knowing the topology of the 
native conformation Thus, the topology itself also dictates the impact of chemistry 

on the folding process. 

The purpose of the present work is three-fold: 

i) determination of the folding bottlenecks for the cellular prion and the identification of 
the key amino acids taking part to the corresponding crucial contacts; 

ii) making connection between the set of such key residues with those that are known to 
be associated with harmful PrP sc mutations. As argued before, a random mutation on key 
positions will usually result in a disruption of the folding process. Only fine-tuned mutations 
can lead to a wild-type-like native state (as for viral enzymes mutating under drug attack 
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|33[ ) or in another viable structure (as postulated by Prusiner for the prion [l|]); 

iii) furthermore, we will focus our attention on how the topological differences between 
the native states of Dpi and PrP c may have impact on the folding process thus aiding or 
avoiding the formation of misfolded conformers. It is important to stress that our study is 
based on the native-state structures of PrP c and Dpi. Hence, although we can confidently 
identify the crucial folding residues, we cannot confirm explicitly that their mutation results 
in a different native state form. 



II. MATERIALS AND METHODS 

A. The model for prion 

Our model builds on a schematic representation of the three-dimensional structure of 
the protein, where the amino acids are replaced by effective centroids identified with the 



C a atoms. The bias towards the native state is introduced through a Go-like []35| energy 
scoring function that rewards the formation of native contacts between the centroids. The 
list of native interactions was compiled from the knowledge of the coordinates of atoms of 
the human prion protein (PDB code lqlx). A pair of amino acids is considered in interaction 
if any pair of their heavy atoms, i and j, have a native separation smaller than the distance, 
1.244(i?j + Rj), where the point of inflection of the Van der Waals interaction occurs (Ri 
denotes the Van der Waals radius of atom i) . The values of Ri are taken from Ref . |3B| . In 
this way one defines a symmetric matrix, known as contact map, A, whose entries, Aij are 
equal to 1 if the i-th and j-th centroids interact, and zero otherwise. 

The energy function for our system includes terms that are routinely used in standard 
molecular dynamics (MD) simulations on biopolymers |37|j38[| . It is composed by a "bonded" 



and a "non-bonded" term: E = Vb + Vnb- The former accounts for the constraints such as 
the peptide bond length and Ramachandran angles bias acting on the aminoacids at a local 
level. Its explicit form is: 



V B = g-V p + h-V a + k-V d , (1) 

N-l 

V P = v £ [r hl+1 - r&i) , (2) 
i=i 

N-2 2 

Va — a ^ (#1,1+1,1+2 — ^"+1,4+2) > (3) 
i=l 
N-3 

V d = aJ2 (l - cos (ji,i + i,i +2 ,i+3 - , (4) 



i=l 



where ry is the distance between residue i and j, 9ij,k is the angle with the j-th amino 
acid as vertex and the i-th and the k-th as edges and r^j^i is the dihedral generated by the 
i-th, the j-th, the k-th and the l-th amino acid. The n as superscript denotes the native 
state value, a represents a suitable scale factor to fix the temperature scale (set to 10 in our 
simulations) which is given in dimensionless units. 

The minimum of expression (|1|) is precisely attained in correspondence of the native 
state. However the cooperative character of the folding process |39| can be captured only 
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by introducing an explicit bias towards the formation of native interactions |[40|| . This is 
accomplished through the second term, Vats, which weights the interaction of any pair of 
non- consecutive centroids with a Van der Waals-like potential: 



N N 



V NB = a £ V lJ = a £ l 5 (^|) - s(lfe) Aij ] . (5) 

i<j-3 i<j-3 V \r i3 - / yr^ / 

This choice rewards native interactions and disfavors non-native ones. The form of the 
energy scoring function ensures that the global minimum is attained in correspondence of 
the native state, regardless of the precise values for the (positive) couplings parameters g, 
h, k. The values used in the present study are: g — 50, h — 5 and k = 0.3; we have checked 



that the results are robust to their variation within a certain range ||34|| . 

The dimensionless unit of time in our MD simulations is 5 • 10~ 3 and conformations 
were sampled every 500 time steps to allow a sufficient uncorrelation at temperatures below 



the folding transition. The multiple-histogram technique f4l| has been applied to reweight 



quantities such as the average internal energy collected in more than 50 equilibrated runs 
at different temperatures in the range 2-8. The optimal reweighting allows to calculate 
thermodynamic quantities, such as the average internal energy and the specific heat for an 
arbitrary temperature, T. 

As visible in Fig. [TjA, the specific heat exhibits a single peak which signals the folding 
transition [l2|. In addition to this overall characterization of the equilibrium thermodynam- 



ics, it is possible to monitor how the individual interactions that are present in the native 
state are formed as the folding progresses. Indeed, although all native contacts are ener- 
getically favored in the same way, their entropic cost of formation may significantly vary 
according to their locality, burial/exposure etc. . 

The crucial contacts for the folding process are identified according to the method given 
in Refs. P3] , |5^ , |i0[] , which singles out the contacts giving the largest contribution to the 
overall specific heat at the folding transition temperature, Tp. 

cij = d ~^ L = {(VnbV^) - (V NB ) {V l3 )) /kT* (6) 

The physical interpretation of this procedure is that these special contacts are precisely 
those that, having a significant and sudden formation at Tp, act as bottlenecks for the folding 
process. Qj measures the sensitivity of the average energy of the system to perturbation 
occurring at the contact (i, j). 

In our studies, a long run carried out at T F of 2 • 10 7 time steps, allowed a detailed 
analysis for all contacts. In addition to the identification of crucial contacts, it is useful to 
introduce a similar measure for the individual amino acids through the single residue specific 
heat (SRSH) resulting from the contribution of all contacts to which a given site, i, takes 
part to: 

Ci = £ ^ij c ij C 7 ) 

j 

SRSH's have been used to identify folding and unfolding phases and to pinpoint residues 
that play a relevant role in the process |3"3"| , |3"1|] . The crucial folding steps can thus be identified 



with two (related) criteria: namely either through the contacts or the sites with the largest 
specific heat of formation. 
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B. Folding pathways 



Another issue that we have addressed in our study is to establish the presence of al- 
ternative pathways taking from unfolded conformations to the native state. Again, this 
is accomplished by examining in detail both the thermodynamic and the kinetic relevant 
features extracted from the model. The comparison of dynamic trajectories connecting un- 
folded and native-like conformations has revealed the presence of several folding events that 
have to occur for the protein to fold. The main ones that emerged from our analysis, and 
that are particularly apt for characterizing the folding process, correspond to the formation 
of contacts involving three sets of contiguous residues. The first set, Bh, from residue 125 
to residue 172 comprises all N-terminal residues (/3-sheet included) before the helix 2 (de- 
fined according to the indications in the PDB file). The second set, hi, from residue 173 to 
194 comprises residues of helix 2, and finally the third set, h2, from residue 200 to 228 is 
composed of residues from helix 3. 

The four observed events are: 

• Event Bh: formation of contacts between residues within Bh, 

• Event Bh-hl: formation of contacts between Bh and hi, 

• Event Bh-h2: formation of contacts between Bh and h2, 

• Event hl-h2: formation of contacts between hi and h2. 

For convenience of presentation, in order to describe the four events on a similar footing, 
we focus on the top 10 contacts (in terms of specific heat contributions c^) in each interacting 
group. It is interesting to note that, among this top contacts, those occurring within an a 
helix almost never occur. This can be explained since local contacts, among others, are easily 
formed, even at high temperature, and hence do not constitute a bottleneck for overcoming 
the configurational barrier to the native state. 

We computed the fraction of formed contacts, Qi (i = Bh, Bh-hl, Bh-h2, hl-h2), within 
each of the 4 groups along the dynamical trajectories of 2 • 10 7 steps at Tp. The natural 
binning size for each of the Qi was 1/10. From this, we obtained a vivid picture of the 
dynamical effects of the free energy landscape at the folding transition (see Fig. |5| and the 
Results section). 

We also carried out 100 ("quenched") folding runs at T = (1 — e)Tp, where e = 1/30, 
slightly below Tp, starting from completely unfolded conformations (thermalised at high 
temperature runs). In all cases, the native state was reached in less than 2 • 10 5 steps. The 
average Q^s as a function of the elapsed time using all 100 trajectories was also calculated 
(Fig. |). 

C. Doppel protein 

We carried out a parallel analysis for the doppel protein, Dpi, whose NMR-resolved 
structure was retrieved from the PDB (PDB code lil7, model 1). For this protein too, we 
built the contact map considering all-atom distances, and obtained both the overall specific 
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heat (Fig. |Tp3 ) and the one associated to individual contacts using the same methods 
described above. A crucial difference in the folding events of the two proteins was revealed 
by monitoring the four events, analogous to those of the previous section. 

For Dpi the sets of amino acids structurally homologous of the prion ones are as follows: 
the set Bh corresponds to residues from 1 to 50 before helix 3; hi to residues from 51 to 
76 (comprising helix 3 and 4) and the set hi to residues from 78 to 96 (comprising helix 5). 
The numeration of the helices is consistent with the secondary structure indications present 
in the PDB file. Free energy profiles and time dependent averages of the reaction coordinate 
are shown in Fig. [| and |5| respectively. 

III. RESULTS 
A. Determination of key residues 

The specific heat curves as a function of the temperature for PrP c and Dpi are shown 
in Figures [l]A and and [l|B, respectively. The single peak in each of them allows to identify 
unambiguously the folding transition in these systems. 

From the long run of 2 • 10 7 time steps performed at Tp, the SRSH's for prion have been 
computed according to eq.(^). We have ranked the residues according to their SRSH (see 
Table 1). It is natural to compare the top key sites from our analysis with those whose 
mutation is known to favor the emergence of PrP sc (highlighted in Table 2). In the top 
10 amino acids of our list we found 4 of the known key mutating sites; whereas in the 
top 30 amino acids their number increases to 12. The statistical relevance of this match 
can be obtained from a combinatorics calculation of the probability to have at least the 
observed number of matches by pure chance. For the top 10 amino acids the probability to 
get a better results by chance is 6.8% whereas for the top 30 amino acids that probability 
becomes 0.02%. Thus, the criterion based on the SRSH allows to identify confidently amino 
acids important for the folding process of the prion protein. 

Their structural location and role is made apparent in Fig. |2] where the prion contact 
map has been colored according to the value of the SRSH. As visible, crucial contacts con- 
nect secondary-structure elements thus forming the protein tertiary structure. In particular 
they connect the (3— sheet-helix- 1 region to the helix-2 and helix-3 [Bh-hl and Bh-h2 sets, 
respectively) and helix-2 - helix-3 (hl-h2 set). The parts of the protein contributing more 
to the specific heat, within our model, are located at the C-terminal strand of the /3-sheet, 
and at the two C-terminal helices, especially helix-3 (see Fig. |3|). 

B. Folding pathways 

The analysis of the long run of 2 • 10 7 time steps carried out at Tp has led to the 
construction of an effective free energy landscape. It is based on the reaction coordinates 
Qi as defined in the Methods section. The free energy as a function of pairs (Qi, Qj) for all 
six combinations of i and j have been calculated as 

F(Q i ,Q j ) = - In (P(Q l ,Q j )) (8) 
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where P{Qi, Qj) is the fraction of times the dynamical trajectory spends in the bin around 
(Qi,Qj)- The contour plots of F are shown in Figures and |4]B for PrP c and Dpi 
respectively. The two energy minima around Q = and in the higher Q region correspond to 
the unfolded and (partially) folded conformations respectively. For prion the plots involving 
set hX-hl show the presence of multiple pathways connecting these two minima. For example, 
the contour plot in the bidimensional space spanned by Qhi-h2 and QBh-hi presents two 
possible successions of events: either the Bh-hl set is formed followed by formation of hl- 
h2 set or vice versa (see Fig. |]A). The same observation can be done for the other plots 
involving set hl-h2. Such pathway ambiguity is absent for the case of doppel (Fig. |]B), 
where trajectories can be grouped essentially in a single set. Indeed, for doppel, on average, 
the formation of the hX-Kl contacts follows the other sets. This result is further confirmed by 
the analysis of the time-dependent average of the reaction coordinates Qi's as obtained from 
the 100 non-equilibrium folding runs. Figures |5)A and p3 show the results for PrP c and Dpi 
respectively. For Dpi the order in the formation of the contacts is the same detected by the 
the equilibrium run: set Bh is formed before set Bh-hl that is followed by Bh-h2 eventually 
followed by hl-h2 interhelical contacts. On the other hand, such definite succession of events 
is not observed for prion. 

IV. DISCUSSION 

The goal of our numerical study was to characterize the folding process of the prion 
protein elucidating the structural mechanisms responsible for the propensity of PrP c to 
mutate into the harmful scrapie form, PrP sc . 

At the heart of our analysis is a topology-based energy function that, by rewarding the 
formation of native-like interactions among the residues, allows to mimic the progressive 
build-up of native structure in a folding process. Remarkably, the progress towards the 
native state does not advance smoothly, but occurs through the overcoming of configurational 
barriers, which are responsible for the appearance of a marked peak in the specific heat curve. 
Crucial contacts (and sites) for the folding process have thus been identified as those giving 
the largest contributions to the specific heat peak (Fig. |2|). 

The top sites isolated by our procedure contain a highly significant fraction of residues 
that are known to enhance the misfolding propensity of PrP c . This indicates that there is 
little freedom in choosing the chemical identity of the amino acids at the sites which have a 
special structural role during the folding process. 

The analysis presented here, based on topological arguments, provides a physically- 
appealing interpretation for the connection between certain site mutations and the arousal 
of the prion disease. Conversely, the key mutating sites not captured by our analysis are 
presumably crucial because they affect some aspects of the process leading to the scrapie 
form that cannot be accounted for by the present method. For example, the alternative 
mechanism leading to misconformation is the template binding of PrP c to PrP sc involving 
the interconversion of the former into the latter one. Consistently with this interpretation, 
amino acids that do not present high SRSH and that are known to be involved in the devel- 
opment of the disease, specifically MET 129, GLU200 and GLU 219, (i.e. mostly charged 
residues with a large exposed surface area). 
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Despite the substantial topological similarities of Dpi and PrP c , the detailed analysis of 
the folding trajectories of both proteins have revealed significant differences of the folding 
pathways. Strikingly, for PrP c the folding pathways can be classified in two main distinct 
routes, whereas for Dpi only one of these is essentially present. For the latter, the series 
of folding steps is more distinctly marked than for the former. Such noticeably different 
behavior can be ascribed to the different set of contacts stabilizing the assembly of the a- 



helices and /3-sheets in the two proteins. It is interesting to note that the contact order GT 



of PrP c is 0.2056 and is lower than that of Dpi which is 0.2298, even though the former has 
somewhat higher number of contacts which is, in turn, reflected by the shift of the folding 
transition peak towards higher temperatures (see Fig. [I]). This difference could explain 
the different average folding time measured in our models for the two proteins. Indeed, 
according to Fig. [|, the sets of contacts in prion are all formed within the first 25000 time 
steps while contacts in doppel contacts can take 4-5 times longer to form. This difference is 
of the same order of magnitude of that predicted by Plaxco on the basis of the contact 
order of the two proteins. Although the present picture would certainly be influenced by 
the introduction of amino acid specific interactions, it is appealing to connect the different 
folding routes to the appearance of misfolded conformers. According to this view, the latter 
would result from following the folding route of PrP c where the inter C-terminal helices 
contacts (set hl-h2) are formed while the contacts with the /3-sheet (set Bh-hl and Bh-h2) 
are not yet formed (Figures |]A and |5)A). Indeed, along this route, which is alternative to 
the one where the two events are interchanged as it occurs in the single route of doppel 
(Fig. (|B and||B), the N-terminal part of the protein is free to rearrange its structure before 
reaching the native-state. This scenario is consistent with several recent models (see e.g. 
Ref. [43,44]) that identify the N-terminal part of the protein as the one that undergoes the 
conformational change that leads the scrapie form of prion. 
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TABLES 



1 

rank 


site 


1 

rank 


site 


1 

rank 


site 


1 


137 


11 


213 


21 


214 


2 


210 


12 


165 


22 


163 


3 


175 


13 


160 


23 


212 


4 


158 


14 


162 


24 


184 


5 


141 


15 


206 


25 


187 


6 


183 


16 


205 


26 


179 


7 


209 


17 


157 


27 


208 


8 


198 


18 


139 


28 


217 


9 


161 


19 


150 


29 


159 


10 


211 


20 


134 


30 


180 



TABLE I. Top 30 sites in prion ranked according to the single residue specific heat (SRSH). 



145, 160, 178, 179, 180, 183, 187, 188, 196, 
198 . 200, 202, 208, 210, 211, 212, 214, 217 



TABLE II. Sites where mutations (in the NMR-solved part of Human PrP) have been observed 
to cause/stop the prion disease or be fundamental for correct folding of PrP c . Sites 214 and 179 
represent cysteine residues that form a disulfide bond needed to observe the PrP c a— form of prion 
[32]. Single residue specific heat (SRSH), as determined by our model, ranking within the top 10 
are in boldface; those ranking within the top 30 are underlined (see Table 1). 
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Figure captions 



1. Figure 1. Temperature dependent specific heat for the prion (A) and doppel (B) 
models. The folding transition temperature is identified at the point where the specific 
heat peak occurs. 

2. Figure 2. Contact map of PrP c (PDB code lqlx) where contacts are highlighted 
according to their specific heat cy of equation (||). Red spots correspond to contacts 
with high Cij, while blue ones correspond low SRSH. The relevant sets of contacts (as 
defined in the text) are framed within the rectangles; (a) set Bh, (b) set Bh-hl, (c) 
set Bh-h2, (d) set hl-h2. Notice that red spots are prevalent in all sets but Bh. 

3. Figure 3. The three-dimensional structure of prion protein where amino acids have 
been colored according to the value of their SRSH. Red spots are located mainly on 
the C-terminal strand of the /3— sheet, on helix-2 and helix-3. 

4. Figure 4. Effective free energy landscape for prion (A) and doppel (B) model. Upper 
and right squares represent the contact maps where the representative set of contacts 
defining the reaction coordinates are highlighted (see Method section). Contour plots 
are shown for each of the six pairs of reaction coordinates (Qi, Qj). White (black) spots 
represent low (high) values of the free energy. The minimum close to (0,0) corresponds 
to the unfolded state whereas the minimum closer to (1,1) corresponds to the folded 
state. In most of the contour plots for prion (A) more than one class of favorable paths 
from the unfolded to the folded conformation is possible whereas for doppel (B) only 
one class of paths is the most favorable. Following these pathways allows to assign a 
preferential order to the events. Thus, on average, the sequence of contact formation 
for doppel is Bh, Bh-hl, Bh-h2 and hl-h2 (see also Fig. |^). For prion, two class of 
equivalent pathways are present: in one case the hl-h2 contacts are formed before the 
other sets of contacts; in the other case h\-h2 contacts follow the formation of all the 
other contacts as for Dpi. 

5. Figure 5. Average time dependence of reactions coordinates Q^s for prion (A) and 
doppel (B) computed on 100 folding runs starting from completely unfolded confor- 
mations. (A) In case of prion the curves corresponding to the different sets (red for 
Bh, green for Bh-hl, blue for Bh-h2 and magenta for hl-h2) present intersections. 
Indeed the average over the different folding pathways results in a non well-defined 
sequence of folding events (B) The folding events in the case of doppel protein have 
a well defined sequence: Bh set (red) is the first to form; then Bh-hl (green); then 
BH-H2 (blue); finally hl-h2 (magenta). 
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FIG. 1. 
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